Cite this Article
Peng Ya, Jiang Zhong-An, Chen Ju-Shi. Surface for methane combustion: O(3P) +CH4 → OH+CH3. Chinese Physics B, 2020, 29(7): 073401  
Permissions
Surface for methane combustion: O(3P) +CH4 → OH+CH3
Peng Ya †, Jiang Zhong-An ‡, Chen Ju-Shi
School of Civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083, China

 

† Corresponding author. E-mail: pengyaustb@sina.com jza1963@263.net

Project supported by the National Natural Science Foundation of China (Grant No. 51574016) and completed while the author was in residence at UNSW, Australia supported by the International Cooperation Training Program for Innovative Talents of USTB.

Abstract

Kinetic investigations including quasi-classical trajectory and canonical unified statistical theory method calculations are carried out on a potential energy surface for the hydrogen-abstraction reaction from methane by atom O(3P). The surface is constructed using a modified Shepard interpolation method. The ab initio calculations are performed at the CCSD(T) level. Taking account of the contribution of inner core electrons to electronic correlation interaction in ab initio electronic structure calculations, modified optimized aug-cc-pCVQZ basis sets are applied to the all-electrons calculations. On this potential energy surface, the triplet oxygen atom attacks methane in a near-collinear H–CH3 direction to form a saddle point with barrier height of 13.55 kcal/mol, which plays a key role in the kinetics of the title reaction. For the temperature range of 298–2500 K, our calculated thermal rate constants for the O(3P) + CH4 → OH + CH3 reaction show good agreement with relevant experimental data. This work provides detailed mechanism of this gas-phase reaction and a theoretical guidance for methane combustion.

PACS: ;34.10.+x;;31.15.xv;;34.50.Lf;
Keyword:;methane combustion;;potential energy surface;;transition state;;kinetic;;ab initio electronic structure calculation;
1. Introduction

Hydrocarbon combustion involves many kinds of chemical reactions, in which reactions of O(3P) with hydrocarbons play an important role in the key initial steps.[1–10] The reaction O(3P) + CH4 → OH + CH3 has been investigated in a lot of experimental and theoretical studies over the decades as a prototype of the hydrocarbon combustion.[11–18] In 2005, all of the available and recommended data was summarized.[19] The thermal rate constants was recommended as k(T) = 7.3×10–19T2.3exp(–3310/T) cm3⋅molecule–1⋅s–1 over the temperature range from 400 K to 2500 K, while the measured rate constants were distributed widely below 400 K. It has also been investigated theoretically by using multifarious approaches.[11,20–23] The Jahn–Teller conical intersection has advanced nature and potential effectiveness.[1,20–22,24–26] Its appearance as the oxygen atom gets close to the hydrogen atom along the H–CH3 bond has not only aroused great theoretical concern but also brought some uncertainty to kinetics. The 3-fold symmetry leads to a 3E state, which could be split into one 3A′ and one 3A″ after the break of C3v symmetry surface.[27–29] Recently our group reported a series of studies to obtain a better understanding of the reaction mechanisms and clarify the role of saddle point. The barrier heights of 3A′ and 3A″ first-order saddle points formed by O(3P) attacking hydrogen in the near-collinear H–CH3 direction are 14.17 and 14.13 kcal/mol, respectively.[4–6] Our calculated imaginary frequency for the TS is 2376i in cm–1, corresponding to a narrow potential barrier.

To achieve accurate kinetic investigation of the title reaction, an accurate potential energy surface (PES) is basically required.[25,26,30–32] A series of ab initio PESs for the title reaction have been constructed. Despite the first full dimensional analytical one (labeled as PES-1998) was constructed by Corchado and co-workers,[24] the most widely used one was constructed by Espinosa-Garcia and co-workers.[33] The PES-1998 was on the basis of studies by Joseph et al.[34] and Jordan et al.,[35] in which PESs were found to be not completely symmetric in all the potential energy terms pertaining to the four hydrogen atoms of CH4. Later, the PES-2000 and PES-2014 for O(3P) + CH4 were reported,[25,33,36] especially the full-dimensional analytical PES-2014 by González–Lavado and co-authors which has been acted as a benchmark to test the accuracy of approximate theories[31,37–40] and can be applied to more complicated reactions. The PES-2014 predicts a barrier height of 14.0 kcal/mol with collinear O–H–CH3 configuration. The ring polymer molecular dynamics (RPMD) method[41] and transition state theory with the least-action tunneling approximation (LAT) were employed to obtain the thermal rate coefficients of this reaction. The RPMD rate theory gives good estimate of quantum statistics in low temperature range and better one as the temperature goes higher.[26,41]

In this article, we present a ground-state PES for the O(3P) + CH4 → OH + CH3 reaction with the 3A″ transition state (TS) in gas phase. This PES is constructed with a modified Shepard interpolation scheme,[42,43] which has been successfully applied in the PES construction of several polyatomic reactions, including H + SiH4, H + CH4 reaction systems.[44–48] Kinetic calculations obtained from this PES are compared in detail with previous available data. The methodology for the construction of PES are described in Section 2. The properties of the fitted surface are discussed and evaluated in Section 3. Section 4 discusses kinetic calculation results. Finally, a summary of the fitted PES is presented in Section 5.

2. Theoretical methodology
2.1. Electronic structure calculations

Extensive computational studies at different levels are performed in this work to investigate the PES of the O(3P) + CH4 → OH + CH3 reaction. The Gaussian 09 suite of ab initio programs[49] are chosen and used to perform all the ab initio calculations. The geometry coordinates, electronic structures symmetries and units for the reaction system defined in our previous work[4–6] are used throughout, unless otherwise specified. One of the major difficulties of the present work is to find a reasonable balance between ab initio accuracy and computational cost to yield smooth PESs and reaction barriers. A recent study by Truhlar et al. pointed out that the ab initio with all electrons taken into account and larger basis sets would give a fairly accurate description of the reaction barrier heights.[50] Good convergence could not be achieved when methods or basis sets are changed since the existence of PES crossing may make the estimate of saddle point unclear. This effect also appears in such a similar system as O(3P) + HF reaction, in which PES crossings originate from multiple states correlating to the triplet oxygen and intersect many times.[51] In fact, a completely satisfactory set for the global PES could neither be found by using higher level method nor by larger basis sets.

In our present work, the coupled-cluster singles and doubles approach plus quasi-perturbative triples excitations method, CCSD(T), is employed as it is suitable for accurate description of polyatomic molecules.[25,33,52] In order to reduce artificial error in the estimate of barrier height and long range region, the contribution of the inner core electron correlation effects for the energies are taken into account by using large and cost basis sets aug-cc-pCVQZ, which are modified and optimized in the form of (12s6p5d4f1g)/[7s6p5d4f1g][33,50] (the discussion of this ab initio scheme can be found elsewhere[4–6]). The paths of chemical reactions are followed by integrating the intrinsic reaction coordinate (IRC) with the Euler steepest-descent theory. The minimum energy path with a total of 196 calculated points is determined through IRC calculation from –1.9 a.u.1/2Bohr on the reactant side to + 1.4 a.u.1/2Bohr on the product side, and it is shown schematically in Fig. 1 in which various stationary points are indicated.

Fig. 1. O(3P) + CH4(X1A1) reaction path obtained at the all-electron CCSD(T)/optACVQZ level. The energy of the reactants is set as zero.
2.2. Form of the PES

The energies, the energy first order derivatives and the energy second order derivatives at a set of data point geometries are employed to construct the surface with the modified Shepard interpolation method, which is proposed by Collins and has been successfully applied to many systems of more than four atoms for systematically generating ab initio PES.[42,43] A 2nd order Taylor series expansions is employed to represent the form of PES. The resulting continuous surface interpolates all the data points with all their corresponding symmetry equivalents in a 3N-dimensional underlying surface landscape, which must be calculated at all data points scattered throughout the whole configuration space.

In the present work, the six-atom reaction system is represented using a general and straightforward system of coordinates given by linear combinations of inverse interatomic distances. It requires 0.5×[6×(6-1)] inverse interatomic distances Qi to describe the positions of all atoms for the PES construction. Therefore, at a given configuration of the title reaction system and its all symmetry-equivalent copies, the optimal coordinates are defined as a vector Q = {Q1,Q2,Q3,…,Q15}. The final approximation to the surface is represented by the modified Shepard interpolation, which sums over all reference data and their symmetry equivalent configurations. The form of 2nd order local Taylor interpolation expansions Ti(Q) centered on the set of Ndata data point geometries is

where

The final surface as a weighted average of these Taylor expansions is

where g is the group element; g ∈ G means the inclusion of all elements g in a group G, which acts to interchange symmetry equivalent atoms; goi means that the ith configuration is transformed by g. Thus, the final and total symmetry of the reaction system is gained from their total direct product. The expressions for the normalized weight function ωgoi(Q) and Tgoi(Q), which weight the contribution of the Taylor expansion for each data point to the final surface energy function in the corresponding configuration space, have been discussed in detail elsewhere.[53,54]

2.3. Development of the surface

The selection of data points in Eq. (3) for the final Shepard interpolation is determined with the GROW2.2 package,[55] in which a set configurations are chosen as an initial start for generating subsequent coming points. This iterative method was developed by Collins and co-workers.[6,42,43,53] These points should provide an accurate picture that gives the accurate description of important chemical and kinetic regions in the reaction reactant side to the product side. Not only should the stationary points, points on minimum energy path and those close to the path be selected, but also more data points should be taken into account in view of physical considerations, particulary, the analysis of potential contour plots. For all these points, the electronic energies as well as the energy first and second order derivatives in all directions are obtained by ab initio calculations. From the initial points along the IRC, the size of data set will rise gradually by using the above iterative procedure. The iterative development in the construction of the surface is very efficient for polyatomic molecule reactions (the discussion of this iterative procedure can be found elsewhere[6,42,43,53]).

3. Properties of the PES

Our surface (denoted as PES-2019) is a 12-dimensional potential energy surface and constructed utilizing the methods illustrated above for hydrogen-abstraction reaction in the ground state from methane by triplet atom O. According to the quasi-classical trajectory (QCT) simulations proposed by Collins and Jordan,[42,43] the number of the data point sets chosen is several hundreds to thousands. The data set for the PES-2019 contains 9017 symmetry unique geometries. As the initial set, 196 geometries come from the IRC and 370 geometries are chosen carefully to accurately represent the region of interest, especially the important areas near the stationary points and van der Waals (vdW) areas. All the other geometries are determined by the QCT calculations. The central differences are employed to evaluate numerical Hessians to produce second derivative of energies at every small step. The QCT simulations concentrate on the geometries with the ∠OHC angle larger than 100° and the geometries near the region of reactant and saddle point (SP). To produce a sample of geometries for describing kinetic relevant region, more than two thousand of molecular configurations are picked up by randomizing the selection in the integration of trajectories generated from the reactant side. The average absolute error and the root-mean-square error (RMSE) from our surface are obtained to be 0.511 and 0.804 kcal/mol, respectively, which are quite smaller compared to the barrier height.

3.1. Weakly bound Interaction

This surface gives an accurate and clear topology of the regions close to the stationary points and the weakly bound interactions. The geometry parameters of various stationary points obtained from ab initio geometry optimizations and from the PES-2019 are collected in Tables 1 and 2, in which the available computational and experimental data from literature are also collected, while Fig. 2 shows the optimized geometries of the saddle point and vdW wells.

Fig. 2. Geometric parameters (distance in Bohr and angles in degree) of the saddle point (SP) and the intermediate complexes: H3CH…O (WellR), O…H3CH () and H3C⋯HO (WellP).
Table 1.

Properties of reactants, products and saddle points on the PES-2019 for in comparison with the previous reported data (all bond length in Bohr, bond angle and dihedral angle in degree, frequency in wave number and energy in kcal/mol).

.
Table 2.

Properties of complexes at reactant side (WellR, ) and product side (WellP) and previous reported values (bond length in Bohr, bond angle and dihedral angle in degree, energy in kcal/mol and frequency in wave number).

.

The PES-2019 surface is a 12-dimensional hypersurface and not easy to be represented graphically. Several two-dimensional contour cut plots of the surface are made to check the reliability of our interpolation. Typical cuts of the PESs are shown in Figs. 3 and 4. The interaction of the O atom with methane leads to more than one vdW wells.[25,33,36] The van der Waals complexes play important roles in reaction dynamics.[56] The PES-2019 reproduces our previous ab initio calculations[4–6] in the long range interaction region to support a vdW potential well H3CH⋯O (WellR) in the entrance channel in Fig. 3, where the contour is plotted as a function of ROH [1.6 Bohr, 3.0 Bohr] and RCH [5.6 Bohr, 6.9 Bohr]. The present PES properly describes the structure of WellR since it has been identified to be reactive by the IRC calculation to connect to the corresponding saddle point. Moreover, several dozens of points around it have been chosen carefully in the initial set. The second vdW potential well O⋯H3CH () is optimized and confirmed by frequency calculation but no IRC is found to connect with. Although additional (nearly one hundred) nearby points are taken into account in our PES construction, the contour plot of regions around the is still not as smooth as that of WellR. The depth of reactant well WellR obtained from the PES-2019 is –0.53 kcal/mol, which corresponds to a vdW interaction between the triplet O atom and CH4. It is weakly bounded but helpful to orient the reaction toward WellR as the O⋯H–C axis.

Fig. 3. Contour plot (kcal/mol) of the PES-2019 in the entrance channel. The energy of the reactants is set as zero.
Fig. 4. Contour plot (kcal/mol) of the PES-2019 in the region around the transition state as functions of RCH and ROH. The energy of the reactants is set as zero.
3.2. Saddle point

A quasi-collinear O⋯H⋯CH3 transition state appears when the triplet oxygen atom gets close to the one of four hydrogen atoms along a carbon hydrogen bond after the reactant well. The saddle point geometry in the title reaction has been somewhat uncertain for a long time.[57] Three saddle points for the CH4 + O(3P) → OH + CH3 reaction have been identified by more extensive electronic structure determinations in our previous work,[4,5] and it plays a key role in the kinetics of the reaction. Under Cs symmetry, two SP1st are identified in two symmetries 3A″ and 3A′, respectively. They get separated by one second-order SP2nd(3E) that arises on account of the C3h symmetry Jahn–Teller coupling. In the present work, the first two lowest saddle points are optimized and have a length of the breaking carbon hydrogen bonds of 2.643 and 2.649 Bohr, respectively. The barrier height of the lowest SP1st(3A″) is 0.04 kcal/mol lower than that of the second lowest one SP1st(3A′). The lowest SP1st(3A″) with the IRC through it and some points around it are taken into account in the initial set. In the surface construction, the second lowest SP1st(3A′) is also taken into account but without IRC or any other points around it. Since one 3A″ and one 3A′ surface propagate through the corresponding saddle point region, respectively, and correlate with the same ground-state products, it is assumed in the present work that the following rate constants are multiplied by 2 times with the consideration of the flux through the relative lower 3A″ surface. This approach has been employed in some previous reports.[22,36] All discussions in this paper refer to the 3A″ surface unless otherwise stated.

A comprehensive comparison of barrier heights from various reported references and our calculations is represented in Fig. 5. The barrier height obtained from our surface is 13.55 kcal/mol, which is close to the best estimated value 14.08 kcal/mol,[25] value 14.0 kcal/mol from the PES-2014[33] and our previous reported theoretical value 14.13 kcal/mol from CASSCF calculation.[4] The previous icMRCI + Q ab initio calculations[4] give a barrier height (without zero-point-energy) much lower than the present CCSD(T) = FULL/optACVQZ, which may due to the following two reasons. First, the previous icMRCI + Q value for the barrier height was obtained by single point energy calculation without full geometry optimization. Second, the CCSD(T) = FULL/aug-cc-pCVXZ (X = D, T, Q) calculations show a dramatic influence of basis set on the description of correction energies for the barrier height. It is clearly demonstrated that the barrier height has a downward tendency with further basis set extension.

Fig. 5. The comparison of barrier heights from various reported references and our calculations. a: from Ref. [11]; b: from Ref. [61]; c: from Ref. [21]; d: from Ref. [22]; e: from Ref. [17]; f: from Ref. [25]; g: from Ref. [24]; h: from Ref. [36]; i: from Ref. [33]; j: present work at the CCSD(T) = FULL/ACVxZ (x = d, t, q, and 5) level; k: present work at the CCSD(T) = FULL/optACVxZ (x = d, t, and q) level; l: present work PES-2019; m: from Ref. [6].

In comparison with the high-accurate icMRCI+Q calculations,[4] the CCSD(T) calculations in the present work give a lower saddle point, which overestimates the rate constant determined by the TST-based methods for the title reaction. On the other hand, the subtle aspects topology derived from the surface based on the CCSD(T) calculations gives a better description than those points (or curve) obtained from icMRCI+Q. The rate constant calculations from our surface tend to average out the above two effects.

Contour plot (kcal/mol) for the PES-2019 in the region around the lower saddle point as functions of RCH and ROH is presented in Fig. 4. The length of the breaking carbon hydrogen bond of transition state on the PES-2019 is 2.643 Bohr, being 25.6% (0.539 Bohr) longer than that of the CH4 molecule at equilibrium geometry as shown in Table 1. In Fig. 5 of Ref. [33], the energy of the PES-2014 rises much slowly but straightly as the system moves from the entrance/exit channel to the saddle point. In our Fig. 4, however, the energy is not slow or direct toward the saddle region. The different surface shapes of the transition state could be caused not only by differences of the barrier heights discussed above, but also by the differences of the imaginary frequencies. The imaginary frequencies at the transition states for the PES-1998, PES-2000, and PES-2014 are 1507 i,[24] 1549 i,[36] and 1919 i[33] cm–1, respectively. It is clear that the PES-2014 has the thinnest potential energy barrier among them. The harmonic frequencies calculated on the PES-2019 agree well with ab initio ones and give an imaginary frequency of 2315 i cm–1.

Although some noticeable discrepancy can be observed in contour plots in Fig. 5 of Ref. [33] and Fig. 4, as well as in the description of the second vdW well O⋯H3CH, both surfaces are broadly similar in the overall region. Besides slightly different aspect of ab initio calculation, the main difference between the PES-2014 and PES-2019 is also caused by different approaches used in surface constructions. The former is constructed using an analytical polynomial which is basically a valence bond-molecular mechanics, while the latter using interpolation method, which may not be suitable for weakly bound interactions but can be improved by adding more points on interested regions. The contour lines in Fig. 4 are smooth and without artificial oscillations in various regions. Although the smooth procedure has little effect on quantum kinetics, it can make a big difference on the QCT calculations.

4. Kinetic calculations

Some quasi-classical trajectory calculations have been performed on the prior version of this surface to study the vibrational distributions.[6] The prior version was constructed with 6100 symmetry unique geometries and it is converged. The main difference is that more points near potential wells were taken into account in the PES-2019 construction to give a better description of the weakly bound interactions. For the pair-correlated CH3(ν = 0) + OH(ν′ = 0) products, the characteristics of angular distributions were analyzed by the QCT calculations. These results including the angular distributions of the OH + CH3 products indicate that the distinct backward scattering are fully stronger than other directions. The angular distribution of the OH product represents a predominantly sharply backward scattering peak relative to the O(3P) beam direction, where no forward scattering peak or any signal from the sideways scattering is presented. This agrees well with the experimental measurement that about 90% of the total available energy remains as OH co-products rotations.[1,57] The QCT calculations on the PES-2019 in the present work and those[6] on its prior version are very similar and no observable divergence is found in the dynamic calculations, so similar dynamics results were excluded in the present manuscript to avoid redundancy. More dynamics properties have been calculated on the PES-2019 and compared with previous theoretical and experimental information.

4.1. Probability of reaction

To evaluate the reaction probability of O(3P) + CH4 → OH + CH3, the QCT calculations are carried out on the PES-2019. For the following calculations, the initial CH4 in its ground rovibrational state is used for all trajectories. The QCT calculations are started with the fragments separation of 30.0 Bohr for the reactants asymptote O(3P) + CH4. In the same way, the trajectories are stopped at products asymptote OH + CH3 of 30.0 Bohr. The integration time step in the above calculation is set as 0.05 fs. Figure 6 shows the reaction probability (derived from one thousand trajectories) as a function of the Ndata. The reaction probability does not provide much clue about the accuracy of the PES but gives some validation to check the convergence of PES. The reaction probability is multiplied by a factor of 4 during calculations since all the hydrogen atoms in methane contribute equally and are treated as distinguishable. Thus, the values of reaction probability in Fig. 6 already include a factor of 4. When the reference data set is grown to about 5000 points by the classical kinetic calculations, our obtained reaction probability does not go up and down eventually, and converges to 0.12 with very small changes. At 300 K, when Ecol is set to be 30.0 kcal/mol for the title reaction, about 600000 trajectories are obtained with bmax of 5.0 Bohr. The maximum impact parameter b is the largest value of the impact parameter, which is determined by increasing its value until no new reactive trajectory is obtained. To test the reliability of these parameters for the QCT calculations, nearly a same number of trajectories are computed at several temperatures including 100 K, 200 K, and 500 K, and these results are very similar.

Fig. 6. The probability of reaction for O(3P) + CH4(X1A1) → OH(X2π) + CH3() represented as a function of Ndata. The dotted line displays the average value for Ndata = 7012, 7440, 7661, 7923, 8333, 8550 and 9017.
4.2. Thermal rate constant prediction

The traditional transition-state theory (denoted as TST) and canonical variational TS (denoted as CVT) theory calculations have been applied into this reaction, as well as canonical unified statistical (denoted as CUS) theory.[4,6] The spin–orbit (SO) coupling of the triplet oxygen atom and singlet methane molecular are integrated in the TST, CVT and CUS theory calculations using the POLYRATE 2015 program, which is a general polyatomic rate constants code and has been successfully employed to many polyatomic molecule reactions.[58] In the mass-weighted Cartesian coordinates, the reaction path starts from the highest point of transition state and goes gradually to the reactants side or in the opposite direction to products side. The effect of tunneling for rate constant calculations is taken into account by using the small-curvature tunneling (denoted as SCT), large-curvature tunneling (denoted as LCT) and LAT.

The thermal rate constants of hydrogen-abstraction reaction have been calculated at different temperatures ranging from 298 K to 2500 K, which are plotted in Fig. 7. It can be concluded from this figure that the thermal rate constants show typical Arrhenius behavior at temperatures of 298–2500 K for the O(3P) + CH4 → OH + CH3 reaction. The newly obtained rate constants at 300 K, 600 K and 1000 K are 8.91 × 10–18 cm3⋅molecule–1⋅s–1,4.32 × 10–15 cm3⋅molecule–1⋅s–1, and 1.22 × 10–13 cm3⋅molecule–1⋅s–1, respectively, by using the basic canonical variational method. Compared with the previously reported values based on the PES-2014 (1.32 × 10–17 cm3⋅molecule–1⋅s–1 with CUS/LAT or 1.50 × 10–17 cm3⋅molecule–1⋅s–1 with RPMD) at 300 K,[59] our rate constant value agrees better with the experimental data from the reported work[12,60] (5.50 × 10–18 cm3⋅molecule–1⋅s–1 –9.00 × 10–18 cm3⋅molecule–1⋅s–1). This rate constant difference derived from the PES-2014 and PES-2019 is due to the different barrier heights as discussed in Section 3.2. The barrier height of our new surface is 0.05 kcal/mol lower than the benchmark value. As we know, the lower the barrier is, the larger the TST rate constants will be. The CVT and CUS are based on the TST method, but differ in the location of the transition state from TST. The other possible relative maxima and re-crossing are also taken into account in the CUS calculations, which reduce the rate constant values.

Fig. 7. Arrhenius plot for the CH4+O(3P) → CH3+OH reaction. RPMD(PES-2014) and CUS/LAT(PES-2014): RPMD and CUS/LAT calculations with PES-2014 from Ref. [59]; Exp. 1: Experimental values from Ref. [60]; Exp. 2: Experimental values from Ref. [12].

The variational effects are induced into CVT calculations by optimizing the dividing surface to minimize the rate constant. The above reasons could account for that the TST rate constants are much larger than those of CVT and CUS.

Under 300 K, our CUS/LAT results are larger than the data of the previous corresponding work obtained on the PES-2014, while at high temperatures ours are smaller. The RPMD rate theory is believed to be exact in the classical high-temperature limit, e.g., 2500 K, it gives the rate constant with the PES-2014 as 4.05 × 10–11 cm3⋅molecule–1⋅s–1, which is very close to the experimental data 4.20 × 10–11 cm3⋅molecule–1⋅s–1. The CUS/LAT calculations on the PES-2014 and PES-2019 give the value of 2.37 × 10–11 cm3⋅molecule–1⋅s–1 and 3.00 × 10–11 cm3⋅molecule–1⋅s–1, respectively. Notwithstanding at low-temperature 300 K and high-temperature 2500 K the calculated rate constants obtained from the PES-2019 are better than those from the PES-2014 using the same method, we cannot deny that the rate constants of the PES-2019 and PES-2014 at the most temperature are pretty close contest. In particular, the rate constant difference between our CUS/LAT on the PES-2019 and RPMD (or CUS/LAT) on the PES-2014 becomes very small as the temperature goes higher, which shows that all quantum mechanical effects in practice become virtually invisible. On the whole, both the results are in basic agreement with the available experimental work.

The ratio of the rate coefficients of the unsubstituted reaction to the deuterated reaction, CH4/CD4, is obtained to analyze the kinetic isotope effects (KIEs) which depend on the masses of reactants. The KIEs at different temperatures obtained from various PESs are listed in Table 3. It has been pointed out that the KIE is less sensitive to the accuracy of the PES, but a very sensitive test of dynamics approaches.[59] As far as we know, not only are the experimentally determined isotope data not available for comparison but also the various theoretical reported values are wildly different. Our calculated KIEs from CVT/SCT and the previous works show similar trends, but their values do not agree well. Although, in general, the present KIEs values obtained from CUS(PES-2019) are lower than the RPMD/CUS calculation on the PES-2000/PES-2014,[59] all these KIEs decrease rapidly as the temperature goes higher in the whole temperature range. This phenomenon is explained by Zhao and co-workers.[31] The quantum tunneling effect of the H atom increases more rapidly than that of the heavier D atom, which makes the difference between the rates of O + CH4 and O + CD4 larger and larger with decreasing temperature.

Table 3.

H/D kinetic isotope effects (KIEs) for the O(3P) + CH4/CD4 reaction.

.
5. Summary

We have carried out kinetic studies of hydrogen-abstraction reaction of CH4 + O(3P) on a 12-dimensional PES. This PES is based on extensive all electrons CCSD(T) method of ab initio theory together with basis sets of modified optimized aug-cc-pCVQZ. The modified Shepard interpolation is used to construct the interpolated PES, of which the best RMSE is 0.804 kcal/mol and all the contour plots are smooth. The classical trajectories are calculated on this surface to check the PES convergence. The canonical variational and statistical theory for thermal rate constant calculations are also carried out on the PES-2019, and it is encouraging to see better agreement between the computed CUS/LAT rate constants and the available experimental measurements. This PES shows that the triplet oxygen atom attacks methane in a near-collinear H–CH3 direction to form a saddle point with barrier height of 13.55 kcal/mol. The property analysis in various regions shows the PES-2019 surface is physically reasonable.

The construction of PES with the 3A″ TS and the corresponding rate constants discussed in this work can contribute to a full understanding of the kinetics of reaction O(3P) + CH4 → OH + CH3. More work is needed on the other reaction channel with the 3A′ TS, which requires the mapping of PES with the 3A′ TS and characterization of the crossing between the two electronic states PESs.

Reference
[1] Zhang J H Liu K P 2011 Chem. Asian J. 6 3132
[2] Jing F Q Cao J W Liu X J Hu Y F Ma H T Bian W S 2016 Chin. J. Chem. Phys. 29 430
[3] Fu B N Shan X Zhang D H Clary D C 2017 Chem. Soc. Rev. 46 7625
[4] Peng Y Jiang Z A Chen J S 2017 J. Phys. Chem. 121 2209
[5] Peng Y Jiang Z A Chen J S 2018 Chin. Phys. 27 023401
[6] Jiang Z A Peng Y Chen J S Lan G Lin H Y 2018 Chin. Phys. 27 063401
[7] Ault B S 2019 J. Mol. Struct. 1176 47
[8] Guan Y F Yarkony D R 2020 J. Phys. Chem. Lett. 11 1848
[9] Liu R Song H W Yang M H 2019 Chinese J Chem. Phys. 32 46
[10] Troya D 2019 J. Phys. Chem. 123 6911
[11] Walch S P Dunning T H Jr 1980 J. Chem. Phys. 72 3221
[12] Baulch D L Cobos C Cox R A et al. 1992 J. Phys. Chem. Ref. Data 21 411
[13] Suzuki T Hirota E 1993 J. Chem. Phys. 98 2387
[14] Wang M L Li Y M Zhang J Z 2001 J. Phys. Chem. 105 2530
[15] Ausfelder F Kelso H Mckendrick K G 2002 Phys. Chem. Chem. Phys 4 473
[16] Troya D Schatz G C Garton D J Brunsvold A L Minton T K 2004 J. Chem. Phys. 120 731
[17] Troya D Garcia-Molina E 2005 J. Phys. Chem. 109 3015
[18] Zhang J M Lahankar S A Garton D J Minton T K Zhang W Q Yang X M 2011 J. Phys. Chem. 115 10894
[19] Baulch D L Bowman C T Cobos C J et al. 2005 J. Phys. Chem. Ref. Data 34 757
[20] González C Mcdouall J J W Schlegel H B 1990 J. Phys. Chem. 94 7467
[21] González M Hernando J Millán J Sayós R 1999 J. Chem. Phys. 110 7326
[22] Roberto-Neto O Machado F B C Truhlar D G 1999 J. Chem. Phys. 111 10046
[23] Shao K J Fu B N Zhang D H 2015 Chin. J. Chem. Phys. 28 403
[24] Corchado J C Espinosa-Garcia J Roberto-Neto O Chuang Y Y Truhlar D G 1998 J. Phys. Chem. 102 4899
[25] Czakó G Bowman J M 2012 Proc. Natl. Acad. Sci. USA 109 7997
[26] Li Y L Suleimanov Y V Green W H Guo H 2014 J. Phys. Chem. 118 1989
[27] Cederbaum L S Domcke W Köppel H 1978 Chem. Phys. 33 319
[28] Domcke W Mishra S Poluyanov L V 2006 Chem. Phys. 322 405
[29] Opalka D Segado M Poluyanov L V Domcke W 2010 Phys. Rev. 81 042501
[30] Czakó G 2014 J. Chem. Phys. 140 231102
[31] Zhao H L Wang W J Zhao Y 2016 J. Phys. Chem. 120 7589
[32] Bowman J M Czakó G Fu B N 2011 Phys. Chem. Chem. Phys. 13 8094
[33] González-Lavado E Corchado J C Espinosa-Garcia J 2014 J. Chem. Phys. 140 064310
[34] Joseph T Steckler R Truhlar D G 1987 J. Chem. Phys. 87 7036
[35] Jordan M J T Gilbert R G 1995 J. Chem. Phys. 102 5669
[36] Espinosa-Garcia J Garcia-Bernaldez J C 2000 Phys. Chem. Chem. Phys. 2 2345
[37] Espinosa-Garcia J 2014 J. Phys. Chem. 118 3572
[38] González-Lavado E Rangel C Espinosa-Garcia J 2014 Phys. Chem. Chem. Phys. 16 8428
[39] Monge-Palacios M González-Lavado E Espinosa-Garcia J 2014 J. Chem. Phys. 141 094307
[40] Jasper A W Sivaramakrishnan R Klippenstein S J 2019 J. Chem. Phys. 150 114112
[41] Suleimanov Y V Aoiz F J Guo H 2016 J. Phys. Chem. 120 8488
[42] Thompson K C Jordan M J T Collins M A 1998 J. Chem. Phys. 108 8302
[43] Morris M Jordan M J T 2014 J. Chem. Phys. 140 204107
[44] Cao J W Zhang Z J Zhang C F Liu K Wang M H Bian W S 2009 Proc. Natl. Acad. Sci. USA 106 13180
[45] Zhang W Q Zhou Y Wu G R et al. 2010 Proc. Natl. Acad. Sci. USA 107 12782
[46] Zhou Y Fu B N Wang C R Collins M A Zhang D H 2011 J. Chem. Phys. 134 064323
[47] Cao J W Zhang Z J Zhang C F Bian W S Guo Y 2011 J. Chem. Phys. 134 024315
[48] Frankcombe T J Collins M A 2011 Phys. Chem. Chem. Phys. 13 8379
[49] Frisch M J Trucks G W Schlegel H B et al. 2010 Gaussian 09, Wallingford CT https://www.gaussian.com/
[50] Zheng J J Zhao Y Truhlar D G 2009 J. Chem. Theory Comput. 5 808
[51] Gomez-Carrasco S Roncero O Gonzalez-Sanchez L Hernandez M L Alvarino J M Paniagua M Aguado A 2005 J. Chem. Phys. 123 114310
[52] Raghavachari K Trucks G W Pople J A Head-Gordon M 1989 Chem. Phys. Lett. 157 479
[53] Bettens R P A Collins M A 1999 J. Chem. Phys. 111 816
[54] Collins M A 2002 Theor. Chem. Acc. 108 313
[55] Jordan M Thompson K Bettens R et al. GROW, version 2.2, a collection of scripts and programs that allow the user to construct molecular potential energy surfaces for either unimolecular/bimolecular reactions or bound-state systems
[56] Cao J W Li F Y Xia W S Bian W S 2019 Chinese J. Chem. Phys. 32 157
[57] Wang F Y Liu K P 2010 Chem. Sci. 1 126
[58] Zheng J, Zhang S, Lynch B J et al. POLYRATE, version 2015, a computer program for the calculation of chemical reaction rates for polyatomics, see https://comp.chem.umn.edu/polyrate/
[59] González-Lavado E Corchado J C Suleimanov Y V Green W H Espinosa-Garcia J 2014 J. Phys. Chem. 118 3243
[60] Cohen N 1986 Int. J. Chem. Kinet. 18 59
[61] Gordon M S Truhlar D G 1986 J. Am. Chem. Soc. 108 5412